library(SingleCellExperiment)
library(scuttle)
library(ggplot2)
library(patchwork)
library(scater)
library(scran)
library(edgeR)
library(limma)
library(pbapply)
library(iCOBRA)
library(harmonicmeanp)
library(stageR)
library(dplyr)
library(ComplexHeatmap)## here() starts at /Users/jg/Desktop/PhD/DD_project/DD_cases/Covid
pdf("./figures/Venndiagram_DD_DE.pdf",
width = 5,
height = 3,
pointsize = 4)
for (j in seq_along(DD_res)) {
levels_status <- levels(pb_bin_ct$`B cell`$Status_on_day_collection_summary)
par(mfrow=c(2,3))
for(i in 1:5){
pval_df <- data.frame(DGE = DGE_res[[j]][[i]]$PValue,
DD = DD_res[[j]][[i]]$PValue)
rownames(pval_df) <- rownames(DGE_res[[j]][[i]])
pad <- pval_df
pad$DGE <- p.adjust(pad$DGE, method = "BH")
pad$DD <- p.adjust(pad$DD, method = "BH")
cobradata <- COBRAData(pval = pval_df, padj=pad)
cobraperf <- calculate_performance(cobradata, splv = "none",maxsplit = 4)
cobraplot <- prepare_data_for_plot(cobraperf, colorscheme = "Dark2", facetted = TRUE)
plot_overlap(cobraplot, main = paste0(names(DD_res)[j], "\n", levels_status[i+1]),
mar = rep(2,4))
}
}
dev.off()## quartz_off_screen
## 2
stageR_res <- vector(mode="list",length = 5)
names(stageR_res) <- names(DD_res)
for (j in seq_along(DD_res)) {
stageR_res_j <- vector(mode="list",length = 5)
names(stageR_res_j) <- levels(pb_bin_ct$`B cell`$Status_on_day_collection_summary)[2:6]
for(i in 1:5){
pScreen <- rep(NA, nrow(DD_res[[j]][[i]]))
for(h in 1:length(pScreen)){
if(is.na(DD_res[[j]][[i]]$PValue[h])){
pScreen[h] <- DGE_res[[j]][[i]]$PValue[h]
} else {
pScreen[h] <- hmp.stat(c(DD_res[[j]][[i]]$PValue[h],
DGE_res[[j]][[i]]$PValue[h]),
w=NULL)
}
}
names(pScreen) <- rownames(DD_res[[j]][[i]])
#confirmation stage
pConfirmation <- as.matrix(cbind(DD_res[[j]][[i]]$PValue, DGE_res[[j]][[i]]$PValue))
dimnames(pConfirmation) <- list(rownames(DD_res[[j]][[i]]), c("DD","DE"))
# stageWise analysis
stageRObj <- stageR(pScreen=pScreen,
pConfirmation=pConfirmation,
pScreenAdjusted=FALSE)
stageRObj <- stageWiseAdjustment(object=stageRObj,
method="none",
alpha=0.05,
allowNA = TRUE)
res <- getResults(stageRObj)
print(colSums(res)) #stage-wise analysis results
if(colSums(res)[1] == 1){
DD_DE_Sig_j_i <- getAdjustedPValues(stageRObj,
onlySignificantGenes=TRUE,
order=FALSE)
} else{
DD_DE_Sig_j_i <- getAdjustedPValues(stageRObj,
onlySignificantGenes=TRUE,
order=TRUE)
}
if(is.null(DD_DE_Sig_j_i)){
DD_DE_Sig_j_i <- NA
}
stageR_res_j[[i]] <- DD_DE_Sig_j_i
}
stageR_res[[j]] <- stageR_res_j
}## padjScreen DD DE
## 0 0 0
## padjScreen DD DE
## 1 1 0
## padjScreen DD DE
## 5 3 5
## padjScreen DD DE
## 1 1 1
## padjScreen DD DE
## 0 0 0
## padjScreen DD DE
## 0 0 0
## padjScreen DD DE
## 95 33 81
## padjScreen DD DE
## 384 213 326
## padjScreen DD DE
## 180 72 138
## padjScreen DD DE
## 15 10 6
## padjScreen DD DE
## 1 1 1
## padjScreen DD DE
## 71 40 55
## padjScreen DD DE
## 326 154 297
## padjScreen DD DE
## 88 52 74
## padjScreen DD DE
## 23 19 9
## padjScreen DD DE
## 0 0 0
## padjScreen DD DE
## 266 161 229
## padjScreen DD DE
## 2345 1580 2239
## padjScreen DD DE
## 1529 812 1474
## padjScreen DD DE
## 1149 453 1125
## padjScreen DD DE
## 0 0 0
## padjScreen DD DE
## 36 21 21
## padjScreen DD DE
## 232 162 159
## padjScreen DD DE
## 99 30 84
## padjScreen DD DE
## 16 14 3
for(j in 1:5){
for(i in 1:5){
if(all(is.na(stageR_res[[j]][[i]]))){
stagewise <- 0
} else if(length(stageR_res[[j]][[i]])==3){
stagewise <- 1
} else {
stagewise <- sum(stageR_res[[j]][[i]][,1] < 0.05, na.rm=TRUE)
}
sig_DD <- rownames(DD_res[[j]][[i]])[which(p.adjust(DD_res[[j]][[i]]$PValue, method="BH") < 0.05)]
sig_DE <- rownames(DGE_res[[j]][[i]])[which(p.adjust(DGE_res[[j]][[i]]$PValue, method="BH") < 0.05)]
separate <- length(unique(sig_DD, sig_DE))
print(c(stagewise, separate))
}
}## [1] 0 0
## [1] 1 2
## [1] 5 3
## [1] 1 1
## [1] 0 6
## [1] 0 0
## [1] 95 8
## [1] 384 191
## [1] 180 50
## [1] 15 11
## [1] 1 1
## [1] 71 36
## [1] 326 110
## [1] 88 43
## [1] 23 27
## [1] 0 0
## [1] 266 135
## [1] 2345 1417
## [1] 1529 555
## [1] 1149 224
## [1] 0 0
## [1] 36 19
## [1] 232 164
## [1] 99 11
## [1] 16 50
par(mfrow=c(2,3))
for (element in pb_bin_ct) {
bin_counts <- assay(element)
of <- colMeans(sweep(bin_counts, 2, element$ncells, "/"))
logitOf <- log(of/(1-of))
plot(x = jitter(as.numeric(element$Status_on_day_collection_summary)),
y = of,
col = element$Site,
pch = 19,
cex = 0.8,
main = "Offset per Status, colored by Site",
xlab = "Covid status")
}visualize_DGE <- function(pb, genes){
cd <- colData(pb)
cd <- droplevels(cd)
libSize <- pb %>%
counts %>%
colSums %>%
unname
gg_list <- lapply(genes, function(gene){
data <- data.frame(cpm = log2(counts(pb)[gene,]+0.5) - log2(libSize+1)+log2(1e6),
status = cd$Status_on_day_collection_summary,
batch = paste0(cd$Site,"_",cd$sex))
data$cpm[is.infinite(data$cpm)] <- NA
levels(data$batch) <- c("Cambridge_female", "Cambridge_male",
"Ncl_female", "Ncl_male")
data$status_batch <- factor(paste0(data$status, "_", data$batch),
levels = paste0(rep(levels(data$status),each=4),
"_",
rep(levels(data$batch), times=2)))
data$dotposition <- data$status_batch
levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
"0.9" = levels(data$status_batch)[2],
"1.1" = levels(data$status_batch)[3],
"1.3" = levels(data$status_batch)[4],
"1.7" = levels(data$status_batch)[5],
"1.9" = levels(data$status_batch)[6],
"2.1" = levels(data$status_batch)[7],
"2.3" = levels(data$status_batch)[8])
data$dotposition <- as.numeric(as.character(data$dotposition))
gg <- ggplot(data = data, aes(x=status,y=cpm)) +
geom_violin() +
geom_jitter(aes(x=dotposition,
y=cpm,
col=batch),
size=2,
width = 0.05) +
ggtitle(paste0("Log CPM: ",gene)) +
theme_bw() +
theme(plot.title = element_text(size=11))
return(gg)
})
return(gg_list)
}
visualize_DD_LOR <- function(pb, genes){
bin_counts <- assay(pb)
cd <- colData(pb)
cd <- droplevels(cd)
of <- colMeans(sweep(bin_counts, 2, cd$ncells, "/"))
odds_global <- of/(1-of)
gg_list <- lapply(genes, function(gene){
odds_gene <- bin_counts[gene,]/(cd$ncells - bin_counts[gene,])
data <- data.frame(odds_gene = odds_gene,
odds_global = odds_global,
LOR = log(odds_gene/odds_global),
status = cd$Status_on_day_collection_summary,
batch = paste0(cd$Site,"_",cd$sex))
data$LOR[is.infinite(data$LOR)] <- NA
levels(data$batch) <- c("Cambridge_female", "Cambridge_male",
"Ncl_female", "Ncl_male")
data$status_batch <- factor(paste0(data$status, "_", data$batch),
levels = paste0(rep(levels(data$status),each=4),
"_",
rep(levels(data$batch), times=2)))
data$dotposition <- data$status_batch
levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
"0.9" = levels(data$status_batch)[2],
"1.1" = levels(data$status_batch)[3],
"1.3" = levels(data$status_batch)[4],
"1.7" = levels(data$status_batch)[5],
"1.9" = levels(data$status_batch)[6],
"2.1" = levels(data$status_batch)[7],
"2.3" = levels(data$status_batch)[8])
data$dotposition <- as.numeric(as.character(data$dotposition))
gg <- ggplot(data = data,aes(x=status, y= LOR)) +
geom_violin() +
geom_jitter(aes(x=dotposition,
y=LOR,
col=batch),
size=2,
width = 0.05) +
ggtitle(paste0("LOR: ", gene)) +
theme_bw() +
theme(plot.title = element_text(size=11))
return(gg)
})
return(gg_list)
}
visualize_DD_logprop <- function(pb, genes){
bin_counts <- assay(pb)
cd <- colData(pb)
cd <- droplevels(cd)
logprop <- log(sweep(bin_counts, 2, cd$ncells, "/"))
gg_list <- lapply(genes, function(gene){
data <- data.frame(logprop = logprop[gene,],
status = cd$Status_on_day_collection_summary,
batch = paste0(cd$Site,"_",cd$sex))
data$logprop[is.infinite(data$logprop)] <- NA
levels(data$batch) <- c("Cambridge_female", "Cambridge_male",
"Ncl_female", "Ncl_male")
data$status_batch <- factor(paste0(data$status, "_", data$batch),
levels = paste0(rep(levels(data$status),each=4),
"_",
rep(levels(data$batch), times=2)))
data$dotposition <- data$status_batch
levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
"0.9" = levels(data$status_batch)[2],
"1.1" = levels(data$status_batch)[3],
"1.3" = levels(data$status_batch)[4],
"1.7" = levels(data$status_batch)[5],
"1.9" = levels(data$status_batch)[6],
"2.1" = levels(data$status_batch)[7],
"2.3" = levels(data$status_batch)[8])
data$dotposition <- as.numeric(as.character(data$dotposition))
gg <- ggplot(data = data,aes(x=status, y= logprop)) +
geom_violin() +
geom_jitter(aes(x=dotposition,
y=logprop,
col=batch),
size=2,
width = 0.05) +
ggtitle(paste0("logprop: ", gene)) +
theme_bw() +
theme(plot.title = element_text(size=11))
return(gg)
})
return(gg_list)
}
visualize_DD_prop <- function(pb, genes){
bin_counts <- assay(pb)
cd <- colData(pb)
cd <- droplevels(cd)
proportions <- sweep(bin_counts, 2, cd$ncells, "/")
gg_list <- lapply(genes, function(gene){
data <- data.frame(proportion = proportions[gene,],
status = cd$Status_on_day_collection_summary,
batch = paste0(cd$Site,"_",cd$sex))
data$proportion[is.infinite(data$proportion)] <- NA
levels(data$batch) <- c("Cambridge_female", "Cambridge_male",
"Ncl_female", "Ncl_male")
data$status_batch <- factor(paste0(data$status, "_", data$batch),
levels = paste0(rep(levels(data$status),each=4),
"_",
rep(levels(data$batch), times=2)))
data$dotposition <- data$status_batch
levels(data$dotposition) <- list("0.7" = levels(data$status_batch)[1],
"0.9" = levels(data$status_batch)[2],
"1.1" = levels(data$status_batch)[3],
"1.3" = levels(data$status_batch)[4],
"1.7" = levels(data$status_batch)[5],
"1.9" = levels(data$status_batch)[6],
"2.1" = levels(data$status_batch)[7],
"2.3" = levels(data$status_batch)[8])
data$dotposition <- as.numeric(as.character(data$dotposition))
gg <- ggplot(data = data,aes(x=status, y= proportion)) +
geom_violin() +
geom_jitter(aes(x=dotposition,
y=proportion,
col=batch),
size=2,
width = 0.05) +
ylim(-0.01,1.01) +
ggtitle(paste0("prop: ", gene)) +
theme_bw() +
theme(plot.title = element_text(size=11))
return(gg)
})
return(gg_list)
}pb_counts_list <- list(Naive_moderate = pb_ct$`naive B cell`[,which(pb_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
Naive_critical = pb_ct$`naive B cell`[,which(pb_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))],
Unswitched_moderate = pb_ct$`unswitched memory B cell`[,which(pb_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
Unswitched_critical = pb_ct$`unswitched memory B cell`[,which(pb_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))])pb_bin_list <- list(Naive_moderate = pb_bin_ct$`naive B cell`[,which(pb_bin_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
Naive_critical = pb_bin_ct$`naive B cell`[,which(pb_bin_ct$`naive B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))],
Unswitched_moderate = pb_bin_ct$`unswitched memory B cell`[,which(pb_bin_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Moderate"))],
Unswitched_critical = pb_bin_ct$`unswitched memory B cell`[,which(pb_bin_ct$`unswitched memory B cell`$Status_on_day_collection_summary %in% c("Healthy", "Critical"))])## logFC logCPM F PValue
## TRBC2 -0.4761722 8.770894 62.98263 1.584941e-12
## ICAM3 -0.5147748 8.565486 60.72267 3.335492e-12
## TGFBR2 -0.5887690 7.925526 59.98712 4.258692e-12
## RNF141 -0.6553944 7.463683 59.21566 5.509029e-12
## MKI67 4.0945860 4.318984 59.42586 5.640242e-12
## SLC43A2 -0.9289432 5.793013 58.13510 7.916771e-12
## MYLIP -0.7271692 7.722048 55.98090 2.607019e-11
## ADD1 -0.4508512 7.924785 52.87000 4.796916e-11
## BLCAP -0.6159960 7.496801 51.76362 7.057862e-11
## RRM2 4.6558042 4.262710 51.41045 1.163980e-10
## logFC logCPM F PValue
## ATP5F1A -0.3454894 8.359434 28.74228 4.368316e-07
## IGHV6-1 1.1558553 4.874379 18.74248 3.229214e-05
## ACTG1 -0.3194227 9.447387 18.32626 4.351886e-05
## GDI2 -0.2332617 8.832662 17.94472 4.629560e-05
## HSPA8 -0.2259126 9.136426 15.85318 1.208046e-04
## HMGB2 -0.3538589 7.858732 14.95412 1.925820e-04
## LITAF -0.2540148 8.096516 14.50821 2.264128e-04
## IGLV1-51 0.9283355 6.242311 14.45082 2.588274e-04
## SHMT2 -0.2586742 8.166360 13.93754 2.964212e-04
## PKM -0.2190923 8.514743 12.81372 5.065802e-04
pdf("./figures/figure3.pdf",
width = 8,
height = 5,
pointsize = 4)
DD_prop_plots[[3]] + DGE_plots[[3]] + plot_layout(guides = "collect")
dev.off()## quartz_off_screen
## 2
## logFC logCPM F PValue
## LINC01238 -1.2922121 2.983519 7.841155 5.999375e-03
## SYNC -0.4995439 5.612165 7.935728 5.713336e-03
## PIK3R5 -0.4212666 5.776195 9.260385 2.906765e-03
## HNRNPUL2 -0.3752241 5.923456 7.463781 7.297005e-03
## HMGB2 -0.3538589 7.858732 14.954121 1.925820e-04
## ATP5F1A -0.3454894 8.359434 28.742279 4.368316e-07
## SAMM50 -0.3201015 6.653094 8.159834 5.090509e-03
## ACTG1 -0.3194227 9.447387 18.326257 4.351886e-05
## DDT -0.2913479 8.224187 10.601145 1.551812e-03
## UBA1 -0.2710409 7.063174 7.961008 5.639294e-03
## logFC logCPM F PValue
## MKI67 4.0121620 4.318984 44.40367 1.050734e-09
## CDKN2D -0.6837345 7.275963 33.51255 6.339351e-08
## PSAT1 3.3397088 3.811498 31.77674 1.269090e-07
## TRBC2 -0.4243709 8.770894 31.51889 1.408042e-07
## RRM2 4.2225069 4.262710 30.89135 2.146758e-07
## BIRC5 3.8207565 4.105887 30.52895 2.319837e-07
## UBE2C 4.3883574 3.738962 27.41210 7.581503e-07
## SEMA4B -0.6655763 7.226843 26.10699 1.367672e-06
## PCLAF 2.9412461 4.277981 25.87200 1.562934e-06
## BUB1 2.6976148 4.102725 25.68458 1.565488e-06
## logFC logCPM F PValue
## IFITM1 0.7048667 9.359272 15.57898 0.0001553255
## IGHV6-1 1.1491385 4.874379 13.20591 0.0004198879
## MCM5 -0.3864310 7.280573 12.76297 0.0005190805
## SNX6 0.2209095 8.361851 12.55092 0.0005748248
## MYC 0.4780873 8.027438 11.92146 0.0008381220
## S100A8 1.1834574 8.245194 11.59541 0.0009842203
## IGKV1D-8 1.1762445 4.358691 11.28444 0.0010632792
## logFC logCPM F PValue
## MCM5 -0.3864310 7.280573 12.76297 0.0005190805
## SNX6 0.2209095 8.361851 12.55092 0.0005748248
## MYC 0.4780873 8.027438 11.92146 0.0008381220
## IFITM1 0.7048667 9.359272 15.57898 0.0001553255
## IGHV6-1 1.1491385 4.874379 13.20591 0.0004198879
## IGKV1D-8 1.1762445 4.358691 11.28444 0.0010632792
## S100A8 1.1834574 8.245194 11.59541 0.0009842203
## logFC logCPM F PValue
## PSME2 0.5507604 9.476988 23.54312 5.383867e-06
## IGKC 0.6453526 9.328480 22.97122 9.375926e-06
## PMAIP1 1.1439169 8.479374 21.65406 2.421569e-05
## KIAA0040 -0.7636898 8.622081 19.16872 3.364242e-05
## NFKBID 1.0154309 8.356222 19.95561 3.416495e-05
## DDX21 0.4502855 9.322344 17.97053 5.154182e-05
## EIF2AK2 1.2710702 8.693864 19.70351 7.317991e-05
## CD27 -0.7000375 8.614283 17.48454 7.418096e-05
## SMIM14 -0.4986694 9.076394 16.29166 1.090757e-04
## ARHGAP30 -0.5941593 8.723185 16.06569 1.207913e-04
## logFC logCPM F PValue
## GPM6A -1.5599708 7.348838 13.811090 0.0004432711
## C3orf38 -1.1731767 7.454788 12.955172 0.0005066061
## MID1IP1 -1.1504175 7.432318 10.821879 0.0014068824
## OSBPL10 -0.9746514 7.706646 13.537482 0.0003856664
## BCL7A -0.8985348 7.871920 12.755253 0.0005953312
## TP53I11 -0.8918054 7.622110 9.614723 0.0025329619
## GPAA1 -0.8884655 7.861906 14.407231 0.0002576068
## RNF113A -0.8447010 7.695519 9.766845 0.0023496663
## TSPYL1 -0.8150866 8.014314 14.454007 0.0002521083
## ALDH16A1 -0.8027842 7.974603 11.063781 0.0013002066
## logFC logCPM F PValue
## EIF5A 0.6137839 9.754504 32.12061 1.512792e-07
## DDX21 0.6480403 9.322344 27.12606 1.083771e-06
## IGKC 0.8428443 9.328480 28.18164 1.322635e-06
## COX6A1 0.4530417 9.825657 21.02758 1.363526e-05
## CD27 -0.9792277 8.614283 19.43624 3.235216e-05
## RHOA 0.4199441 10.007795 18.86143 3.482467e-05
## CD1C -1.2118898 8.643620 20.93737 3.570295e-05
## ZNF581 -1.1502640 8.320019 18.66492 3.795786e-05
## TNIP2 1.2475800 7.881964 18.31475 4.427646e-05
## PSME2 0.5870198 9.476988 18.50603 4.451462e-05
## logFC logCPM F PValue
## LCN8 -5.9292313 6.824761 13.91315 6.293238e-04
## GPM6A -2.2523892 7.348838 13.05752 6.154990e-04
## FCGRT -1.2570043 7.717782 13.08554 4.765082e-04
## CD1C -1.2118898 8.643620 20.93737 3.570295e-05
## ZNF581 -1.1502640 8.320019 18.66492 3.795786e-05
## CD27 -0.9792277 8.614283 19.43624 3.235216e-05
## RIN3 -0.7505344 8.341573 12.68901 5.742842e-04
## FXR1 -0.7379673 8.517854 12.39759 6.591417e-04
## FCRL2 -0.7276242 8.848491 15.26450 1.795342e-04
## LINC01857 -0.6606280 9.070504 13.87236 3.627673e-04
sce_sub <- sce[,sce$cell_type_curated == "naive B cell"]
sce_sub <- sce_sub[,sce_sub$donor_id %in% pb_bin_list$Naive_moderate$donor_id]
lnc <- logNormCounts(sce_sub,name="logNormCounts")
lnc$Site_sex <- paste0(lnc$Site, "_", lnc$sex)gene <- "ACTG1"
DD_prop_plots <- visualize_DD_prop(pb = pb_bin_list$Naive_moderate,
genes = gene)
DGE_plots <- visualize_DGE(pb = pb_counts_list$Naive_moderate,
genes = gene)
print(DD_prop_plots[[1]] + DGE_plots[[1]] + plot_layout(guides = "collect"))df_hlp <- data.frame(counts = assays(sce_sub)$counts[gene,],
lnc = assays(lnc)$logNormCounts[gene,],
Site_sex = lnc$Site_sex,
Status = lnc$Status_on_day_collection_summary,
Donor = as.factor(as.character(lnc$donor_id)))
df_hlp <- df_hlp[df_hlp$Site_sex == "Ncl_male",]
gg1 <- ggplot(data = df_hlp, aes(x=Donor,y=counts, fill=Status)) +
geom_violin() +
facet_grid(~Status, space = "free_x", scales = "free_x") +
ggtitle(paste0("Counts: ",gene)) +
theme_bw() +
theme(plot.title = element_text(size=11),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
gg2 <- ggplot(data = df_hlp, aes(x=Donor,y=lnc, fill=Status)) +
geom_violin() +
facet_grid(~Status, space = "free_x", scales = "free_x") +
ggtitle(paste0("LogNormCounts: ",gene)) +
theme_bw() +
theme(plot.title = element_text(size=11),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))## quartz_off_screen
## 2
# B cell
sce_hlp <- sce[,sce$cell_type_curated == "B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`B cell`$donor_id]
a_entries <- nrow(sce_hlp)*ncol(sce_hlp)
a_0 <- sum(assay(sce_hlp) == 0)
a_1 <- sum(assay(sce_hlp) == 1)
a_2 <- sum(assay(sce_hlp) == 2)
a_r <- sum(assay(sce_hlp) > 2)
# class switched memory B cell
sce_hlp <- sce[,sce$cell_type_curated == "class switched memory B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`class switched memory B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`class switched memory B cell`$donor_id]
b_entries <- nrow(sce_hlp)*ncol(sce_hlp)
b_0 <- sum(assay(sce_hlp) == 0)
b_1 <- sum(assay(sce_hlp) == 1)
b_2 <- sum(assay(sce_hlp) == 2)
b_r <- sum(assay(sce_hlp) > 2)
# immature B cell
sce_hlp <- sce[,sce$cell_type_curated == "immature B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`immature B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`immature B cell`$donor_id]
c_entries <- nrow(sce_hlp)*ncol(sce_hlp)
c_0 <- sum(assay(sce_hlp) == 0)
c_1 <- sum(assay(sce_hlp) == 1)
c_2 <- sum(assay(sce_hlp) == 2)
c_r <- sum(assay(sce_hlp) > 2)
# naive B cell
sce_hlp <- sce[,sce$cell_type_curated == "naive B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`naive B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`naive B cell`$donor_id]
d_entries <- nrow(sce_hlp)*ncol(sce_hlp)
d_0 <- sum(assay(sce_hlp) == 0)
d_1 <- sum(assay(sce_hlp) == 1)
d_2 <- sum(assay(sce_hlp) == 2)
d_r <- sum(assay(sce_hlp) > 2)
# unswitched memory B cell
sce_hlp <- sce[,sce$cell_type_curated == "unswitched memory B cell"]
sce_hlp <- sce_hlp[rownames(pb_ct$`unswitched memory B cell`),]
sce_hlp <- sce_hlp[,sce_hlp$donor_id %in% pb_ct$`unswitched memory B cell`$donor_id]
e_entries <- nrow(sce_hlp)*ncol(sce_hlp)
e_0 <- sum(assay(sce_hlp) == 0)
e_1 <- sum(assay(sce_hlp) == 1)
e_2 <- sum(assay(sce_hlp) == 2)
e_r <- sum(assay(sce_hlp) > 2)## [1] 86.09558
## [1] 9.386102
## [1] 1.924958
## [1] 2.593364
## R version 4.2.0 alpha (2022-04-04 r82084)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ComplexHeatmap_2.14.0 dplyr_1.1.2
## [3] stageR_1.20.0 harmonicmeanp_3.0
## [5] FMStable_0.1-4 iCOBRA_1.26.0
## [7] pbapply_1.7-2 edgeR_3.40.2
## [9] limma_3.54.2 scran_1.26.2
## [11] scater_1.26.1 patchwork_1.1.2
## [13] ggplot2_3.4.2 scuttle_1.8.4
## [15] SingleCellExperiment_1.20.1 SummarizedExperiment_1.28.0
## [17] Biobase_2.58.0 GenomicRanges_1.50.2
## [19] GenomeInfoDb_1.34.9 IRanges_2.32.0
## [21] S4Vectors_0.36.2 BiocGenerics_0.44.0
## [23] MatrixGenerics_1.10.0 matrixStats_1.0.0
##
## loaded via a namespace (and not attached):
## [1] ggbeeswarm_0.7.2 colorspace_2.1-0
## [3] rjson_0.2.21 ellipsis_0.3.2
## [5] rprojroot_2.0.3 circlize_0.4.15
## [7] bluster_1.8.0 XVector_0.38.0
## [9] GlobalOptions_0.1.2 BiocNeighbors_1.16.0
## [11] clue_0.3-64 rstudioapi_0.14
## [13] farver_2.1.1 ggrepel_0.9.3
## [15] DT_0.28 fansi_1.0.4
## [17] codetools_0.2-19 sparseMatrixStats_1.10.0
## [19] doParallel_1.0.17 cachem_1.0.8
## [21] knitr_1.43 jsonlite_1.8.7
## [23] cluster_2.1.4 png_0.1-8
## [25] shinydashboard_0.7.2 shiny_1.7.4.1
## [27] compiler_4.2.0 dqrng_0.3.0
## [29] Matrix_1.5-4.1 fastmap_1.1.1
## [31] cli_3.6.1 later_1.3.1
## [33] BiocSingular_1.14.0 htmltools_0.5.5
## [35] tools_4.2.0 rsvd_1.0.5
## [37] igraph_1.5.0.1 gtable_0.3.3
## [39] glue_1.6.2 GenomeInfoDbData_1.2.9
## [41] reshape2_1.4.4 Rcpp_1.0.11
## [43] jquerylib_0.1.4 vctrs_0.6.3
## [45] iterators_1.0.14 DelayedMatrixStats_1.20.0
## [47] xfun_0.39 stringr_1.5.0
## [49] beachmat_2.14.2 mime_0.12
## [51] lifecycle_1.0.3 irlba_2.3.5.1
## [53] statmod_1.5.0 zlibbioc_1.44.0
## [55] scales_1.2.1 promises_1.2.0.1
## [57] shinyBS_0.61.1 parallel_4.2.0
## [59] RColorBrewer_1.1-3 yaml_2.3.7
## [61] gridExtra_2.3 UpSetR_1.4.0
## [63] sass_0.4.7 stringi_1.7.12
## [65] highr_0.10 foreach_1.5.2
## [67] ScaledMatrix_1.6.0 BiocParallel_1.32.6
## [69] shape_1.4.6 rlang_1.1.1
## [71] pkgconfig_2.0.3 bitops_1.0-7
## [73] evaluate_0.21 lattice_0.21-8
## [75] ROCR_1.0-11 labeling_0.4.2
## [77] htmlwidgets_1.6.2 tidyselect_1.2.0
## [79] here_1.0.1 plyr_1.8.8
## [81] magrittr_2.0.3 R6_2.5.1
## [83] generics_0.1.3 metapod_1.6.0
## [85] DelayedArray_0.24.0 DBI_1.1.3
## [87] pillar_1.9.0 withr_2.5.0
## [89] RCurl_1.98-1.12 tibble_3.2.1
## [91] crayon_1.5.2 utf8_1.2.3
## [93] rmarkdown_2.23 viridis_0.6.4
## [95] GetoptLong_1.0.5 locfit_1.5-9.8
## [97] digest_0.6.33 xtable_1.8-4
## [99] httpuv_1.6.11 munsell_0.5.0
## [101] beeswarm_0.4.0 viridisLite_0.4.2
## [103] vipor_0.4.5 bslib_0.5.0